Practical - Week 3
2024-10-14
We will explore the potential drivers of species distribution.
Exercises:
WorldClim Is a database of high spatial resolution global weather and climate data. You can download gridded weather and climate data for historical (near current) and future conditions.
CHELSA (Climatologies at high resolution for the earth’s land surface areas) is a very high resolution (30 arc sec, ~1km) global downscaled climate data set. It includes climate layers for various time periods and variables, ranging from the Last Glacial Maximum, to the present, to several future scenarios.
A global, spatially explicit characterization of 47 terrestrial habitat types, as defined in the International Union for Conservation of Nature (IUCN) habitat classification scheme.
Habitatmapping: https://github.com/Martin-Jung/Habitatmapping
The RESOLVE Ecoregions dataset, updated in 2017, offers a depiction of the 846 terrestrial ecoregions that represent our living planet.
The CORINE Land Cover (CLC) consists of an inventory of land cover in 44 classes (from Europe). The inventory was initiated in 1985. Updates have been produced in 2000, 2006, 2012, and 2018.
Terra MODIS and Aqua MODIS satellites are viewing the entire Earth’s surface every 1 to 2 days, acquiring data in 36 spectral bands or groups of wavelengths.
MODIStsp: https://github.com/ropensci/MODIStsp/
SEDAC, the Socioeconomic Data and Applications Center, focuses on human interactions in the environment and seeks to develop and operate applications that support the integration of socioeconomic and earth science data.
Steps:
raster)Please download the following file and store it on you data folder:
amphibians_grid.rdsTo create this file, I followed steps from the previous practicals. First I downloaded and cleaned data from GBIF, and then I created a grid cell for Peru of 50 x 50 km.
Check download_amphibiansPE_data_from_GBIF.R to see how the data were downloaded, and preparing_amphibiansPE_grid.R to see how I created the grid
This is how the data look like
We will use sf, rnaturalearth, and tidyverse.
We will also use terra to manage and analyse rasters and vectors, tmap for plotting maps, geodata to get multiple sources of spatial data, and MODIStsp for processing time series from MODIS data.
We will need,
Let’s see how it looks.
To get a planar projection of the area (visit https://projectionwizard.org/).
proj_peru <- 'PROJCS["ProjWiz_Custom_Transverse_Cylindrical_Equal_Area",
GEOGCS["GCS_WGS_1984",
DATUM["D_WGS_1984",
SPHEROID["WGS_1984",6378137.0,298.257223563]],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]],
PROJECTION["Transverse_Cylindrical_Equal_Area"],
PARAMETER["False_Easting",0.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",-74.7949219],
PARAMETER["Scale_Factor",1.0],
PARAMETER["Latitude_Of_Origin",0.0],
UNIT["Meter",1.0]]'We will use it later.
raster)raster)When working with rasters, there are three basic procedures we can recommend you:
extent you need. This is a square cut.projection. But be aware that by reprojecting, you are introducing an additional source of error.NA all the cells outside your polygon(s) of interest.raster)We will create a variable with the spatial extent of Peru to use it to crop the raster data.
Now let’s create an object with the extent.
raster)We will also set the path were we are going to save data.
raster): climateLet’s first download the WorldClim data (version 2.1 climate data for 1970-2000):
bio1)bio12)
You can download these data from https://www.worldclim.org/data/worldclim21.html
or use the package geodata.
raster): climateWe will use the function worldclim_country() from geodata.
class : SpatRaster
dimensions : 2220, 1560, 12 (nrow, ncol, nlyr)
resolution : 0.008333333, 0.008333333 (x, y)
extent : -81.5, -68.5, -18.5, 0 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source : PER_wc2.1_30s_tavg.tif
names : PER_w~avg_1, PER_w~avg_2, PER_w~avg_3, PER_w~avg_4, PER_w~avg_5, PER_w~avg_6, ...
min values : -5.4, -5.4, -5.2, -5.8, -6.9, -8.6, ...
max values : 27.4, 27.3, 27.3, 27.3, 27.0, 26.8, ...
raster): climateNow get the precipitation data from Worldclim.
And filter bio12.
class : SpatRaster
dimensions : 2220, 1560, 1 (nrow, ncol, nlyr)
resolution : 0.008333333, 0.008333333 (x, y)
extent : -81.5, -68.5, -18.5, 0 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source : PER_wc2.1_30s_bio.tif
name : wc2.1_30s_bio_12
min value : 0
max value : 6256
raster): elevationWe will use the function elevation_30s() from geodata. These data come from the Shuttle Radar Topographic Mission and are aggregated to 1km resolution.
class : SpatRaster
dimensions : 2232, 1560, 1 (nrow, ncol, nlyr)
resolution : 0.008333333, 0.008333333 (x, y)
extent : -81.5, -68.5, -18.5, 0.1 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source : PER_elv_msk.tif
name : PER_elv_msk
min value : -25
max value : 6547
raster): land-coverNow, let’s get the landcover data. For this we will use the package MODIStsp.
MODIStsp has not been maintained
The package may not be working properly.
[1] "Surf_Ref_8Days_500m (M*D09A1)"
[2] "Surf_Ref_Daily_005dg (M*D09CMG)"
[3] "Surf_Ref_Daily_500m (M*D09GA)"
[4] "Surf_Ref_Daily_250m (M*D09GQ)"
[5] "Surf_Ref_8Days_250m (M*D09Q1)"
[6] "Snow_Cov_Daily_500m (M*D10A1)"
[7] "Snow_Cov_8-Day_500m (M*D10_A2)"
[8] "Snow_Cov_Day_0.05Deg (M*D10C1)"
[9] "Snow_Cov_8-Day0.05Deg CMG (M*D10C2)"
[10] "Snow_Cov_Month_0.05Deg CMG (M*D10CM)"
[11] "Surf_Temp_Daily_005dg (M*D11C1)"
[12] "Surf_Temp_Daily_1Km (M*D11A1)"
[13] "Surf_Temp_8Days_1Km (M*D11A2)"
[14] "Surf_Temp_Daily_GridSin (M*D11B1)"
[15] "Surf_Temp_8Days_GridSin (M*D11B2)"
[16] "Surf_Temp_Monthly_GridSin (M*D11B3)"
[17] "Surf_Temp_8Days_005dg (M*D11C2)"
[18] "Surf_Temp_Monthly_005dg (M*D11C3)"
[19] "LST_3band_emissivity_Daily_1km (M*D21A1D)"
[20] "LST_3band_emissivity_Daily_1km_night (M*D21A1N)"
[21] "LST_3band_emissivity_8day_1km (M*D21A2)"
[22] "BRDF_Albedo_ModelPar_Daily_500m (MCD43A1)"
[23] "BRDF_Albedo_Quality_Daily_500m (MCD43A2)"
[24] "Albedo_Daily_500m (MCD43A3)"
[25] "BRDF_Adj_Refl_Daily_500m (MCD43A4)"
[26] "BRDF_Albedo_ModelPar_Daily_005dg (MCD43C1)"
[27] "BRDF_Albedo_Quality_Daily_005dg (MCD43C2)"
[28] "Albedo_Daily_005dg (MCD43C3)"
[29] "BRDF_Adj_Refl_16Day_005dg (MCD43C4)"
[30] "AlbPar_1_B1_Daily_30ArcSec (MCD43D01)"
[31] "AlbPar_2_B1_Daily_30ArcSec (MCD43D02)"
[32] "AlbPar_3_B1_Daily_30ArcSec (MCD43D03)"
[33] "AlbPar_1_B2_Daily_30ArcSec (MCD43D04)"
[34] "AlbPar_2_B2_Daily_30ArcSec (MCD43D05)"
[35] "AlbPar_3_B2_Daily_30ArcSec (MCD43D06)"
[36] "AlbPar_1_B3_Daily_30ArcSec (MCD43D07)"
[37] "AlbPar_2_B3_Daily_30ArcSec (MCD43D08)"
[38] "AlbPar_3_B3_Daily_30ArcSec (MCD43D09)"
[39] "AlbPar_1_B4_Daily_30ArcSec (MCD43D10)"
[40] "AlbPar_2_B4_Daily_30ArcSec (MCD43D11)"
[41] "AlbPar_3_B4_Daily_30ArcSec (MCD43D12)"
[42] "AlbPar_1_B4_Daily_30ArcSec (MCD43D13)"
[43] "AlbPar_2_B4_Daily_30ArcSec (MCD43D14)"
[44] "AlbPar_3_B4_Daily_30ArcSec (MCD43D15)"
[45] "AlbPar_1_B5_Daily_30ArcSec (MCD43D16)"
[46] "AlbPar_2_B5_Daily_30ArcSec (MCD43D17)"
[47] "AlbPar_3_B5_Daily_30ArcSec (MCD43D18)"
[48] "AlbPar_1_B6_Daily_30ArcSec (MCD43D19)"
[49] "AlbPar_2_B6_Daily_30ArcSec (MCD43D20)"
[50] "AlbPar_3_B6_Daily_30ArcSec (MCD43D21)"
[51] "AlbPar_1_Vis_Daily_30ArcSec (MCD43D22)"
[52] "AlbPar_2_Vis_Daily_30ArcSec (MCD43D23)"
[53] "AlbPar_3_Vis_Daily_30ArcSec (MCD43D24)"
[54] "AlbPar_1_NIR_Daily_30ArcSec (MCD43D25)"
[55] "AlbPar_2_NIR_Daily_30ArcSec (MCD43D26)"
[56] "AlbPar_3_NIR_Daily_30ArcSec (MCD43D27)"
[57] "AlbPar_1_SWIR_Daily_30ArcSec (MCD43D28)"
[58] "AlbPar_2_SWIR_Daily_30ArcSec (MCD43D29)"
[59] "AlbPar_3_SWIR_Daily_30ArcSec (MCD43D30)"
[60] "BRDF_Albedo_Quality_Daily_30ArcSec (MCD43D31)"
[61] "BRDF_Albedo_SolNoon_Daily_30ArcSec (MCD43D32)"
[62] "Alb_ValObs_B1_Daily_30ArcSec (MCD43D33)"
[63] "Alb_ValObs_B2_Daily_30ArcSec (MCD43D34)"
[64] "Alb_ValObs_B3_Daily_30ArcSec (MCD43D35)"
[65] "Alb_ValObs_B4_Daily_30ArcSec (MCD43D36)"
[66] "Alb_ValObs_B5_Daily_30ArcSec (MCD43D37)"
[67] "Alb_ValObs_B6_Daily_30ArcSec (MCD43D38)"
[68] "Alb_ValObs_B7_Daily_30ArcSec (MCD43D39)"
[69] "BRDF_Albedo_Snow_Daily_30ArcSec (MCD43D40)"
[70] "BRDF_Alb_Unc_Daily_30ArcSec (MCD43D41)"
[71] "BRDF_Alb_BSA_B1_Daily_30ArcSec (MCD43D42)"
[72] "BRDF_Alb_BSA_B2_Daily_30ArcSec (MCD43D43)"
[73] "BRDF_Alb_BSA_B3_Daily_30ArcSec (MCD43D44)"
[74] "BRDF_Alb_BSA_B4_Daily_30ArcSec (MCD43D45)"
[75] "BRDF_Alb_BSA_B5_Daily_30ArcSec (MCD43D46)"
[76] "BRDF_Alb_BSA_B6_Daily_30ArcSec (MCD43D47)"
[77] "BRDF_Alb_BSA_B7_Daily_30ArcSec (MCD43D48)"
[78] "BRDF_Alb_BSA_Vis_Daily_30ArcSec (MCD43D49)"
[79] "BRDF_Alb_BSA_NIR_Daily_30ArcSec (MCD43D50)"
[80] "BRDF_Alb_BSA_SWIR_Daily_30ArcSec (MCD43D51)"
[81] "BRDF_Alb_WSA_B1_Daily_30ArcSec (MCD43D52)"
[82] "BRDF_Alb_WSA_B2_Daily_30ArcSec (MCD43D53)"
[83] "BRDF_Alb_WSA_B3_Daily_30ArcSec (MCD43D54)"
[84] "BRDF_Alb_WSA_B4_Daily_30ArcSec (MCD43D55)"
[85] "BRDF_Alb_WSA_B5_Daily_30ArcSec (MCD43D56)"
[86] "BRDF_Alb_WSA_B6_Daily_30ArcSec (MCD43D57)"
[87] "BRDF_Alb_WSA_B7_Daily_30ArcSec (MCD43D58)"
[88] "BRDF_Alb_WSA_Vis_Daily_30ArcSec (MCD43D59)"
[89] "BRDF_Alb_WSA_NIR_Daily_30ArcSec (MCD43D60)"
[90] "BRDF_Alb_WSA_SWIR_Daily_30ArcSec (MCD43D61)"
[91] "BRDF_Albedo_NBAR_Band1_Daily_30ArcSec (MCD43D62)"
[92] "BRDF_Albedo_NBAR_Band2_Daily_30ArcSec (MCD43D63)"
[93] "BRDF_Albedo_NBAR_Band3_Daily_30ArcSec (MCD43D64)"
[94] "BRDF_Albedo_NBAR_Band4_Daily_30ArcSec (MCD43D65)"
[95] "BRDF_Albedo_NBAR_Band5_Daily_30ArcSec (MCD43D66)"
[96] "BRDF_Albedo_NBAR_Band6_Daily_30ArcSec (MCD43D67)"
[97] "BRDF_Albedo_NBAR_Band7_Daily_30ArcSec (MCD43D68)"
[98] "Vegetation_Indexes_16Days_500m (M*D13A1)"
[99] "Vegetation_Indexes_16Days_1Km (M*D13A2)"
[100] "Vegetation_Indexes_Monthly_1Km (M*D13A3)"
[101] "Vegetation_Indexes_16Days_005dg (M*D13C1)"
[102] "Vegetation_Indexes_Monthly_005dg (M*D13C2)"
[103] "Vegetation Indexes_16Days_250m (M*D13Q1)"
[104] "LAI_8Days_500m (MCD15A2H)"
[105] "LAI_4Days_500m (MCD15A3H)"
[106] "LAI_8Days_500m (M*D15A2H)"
[107] "Net_ET_8Day_500m (M*D16A2)"
[108] "Net_ETgf_8Day_500m (M*D16A2GF)"
[109] "Net_ETgf_Yearly_500m (M*D16A3GF)"
[110] "Gross_PP_8Days_500m (M*D17A2H)"
[111] "Gross_PP_GapFil_8Days_500m (M*D17A2HGF)"
[112] "Net_PP_GapFil_Yearly_500m (M*D17A3HGF)"
[113] "Veg_Cont_Fields_Yearly_250m (MOD44B)"
[114] "Land_Wat_Mask_Yearly_250m (MOD44W)"
[115] "Burned_Monthly_500m (MCD64A1)"
[116] "ThermalAn_Fire_Daily_1Km (M*D14A1)"
[117] "ThermalAn_Fire_8Days_1Km (M*D14A2)"
[118] "LandCover_Type_Yearly_005dg (MCD12C1)"
[119] "LandCover_Type_Yearly_500m (MCD12Q1)"
[120] "LandCover_Dynamics_Yearly_500m (MCD12Q2)"
[121] "Dwnwrd_Srw_Rad_3h_005dg (MCD18A1)"
[122] "Dwnwrd_PAR_3h_005dg (MCD18A2)"
[123] "MAIA_Land_Surf_BRF (MCD19A1)"
[124] "MAIA_Land_AOT_daily (MCD19A2)"
raster): land-coverTo download MODIS data through the NASA http server, we need to create a profile at https://urs.earthdata.nasa.gov/home to get a user and password.
Today, I’ll present you an example and provide you with the data already processed, but you can do this at home.
We will use the Land Cover Type Yearly L3 Global 500m, and download data for Peru from the year 2020.
raster): land-coverWith MODIStsp_get_prodlayers() you can see all the layers of the product:
There are five different land cover classification schemes, we will be using the primary land cover scheme (LC1) which identifies 17 classes defined by the IGBP (International Geosphere-Biosphere Programme), including 11 natural vegetation classes, 3 human-altered classes, and 3 non-vegetated classes.
raster): land-coverWe will use the function MODIStsp() to download data:
MODIStsp(
...,
gui = TRUE,
out_folder = NULL,
out_folder_mod = NULL,
opts_file = NULL,
selprod = NULL,
prod_version = NULL,
bandsel = NULL,
quality_bandsel = NULL,
indexes_bandsel = NULL,
sensor = NULL,
download_server = NULL,
downloader = NULL,
user = NULL,
password = NULL,
download_range = NULL,
start_date = NULL,
end_date = NULL,
spatmeth = NULL,
start_x = NULL,
end_x = NULL,
start_y = NULL,
end_y = NULL,
bbox = NULL,
spafile = NULL,
out_projsel = NULL,
output_proj = NULL,
out_res_sel = NULL,
out_res = NULL,
resampling = NULL,
reprocess = NULL,
delete_hdf = NULL,
nodata_change = NULL,
scale_val = NULL,
ts_format = NULL,
out_format = NULL,
compress = NULL,
test = NULL,
n_retries = 5,
verbose = TRUE,
parallel = TRUE
)raster): land-coverMODIStsp(
gui = FALSE,
out_folder = "data/",
out_folder_mod = "data/",
selprod = "LandCover_Type_Yearly_500m (MCD12Q1)",
prod_version = "061",
bandsel = "LC1",
sensor = "Terra",
user = your_user, # your username for NASA http server
password = your_pass, # your password for NASA http server
start_date = "2020.01.01",
end_date = "2020.12.31",
verbose = TRUE,
bbox = extent, # bbox covering Peru
spatmeth = "bbox",
out_format = "GTiff",
compress = "LZW",
out_projsel = "User Defined",
output_proj = 4326,
delete_hdf = TRUE,
parallel = TRUE
)raster): land-coverLet’s read the layer and process it.
class : SpatRaster
dimensions : 4800, 3360, 1 (nrow, ncol, nlyr)
resolution : 0.004166667, 0.004166667 (x, y)
extent : -82, -68, -19, 1 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source(s) : memory
varname : MCD12Q1.061_LC_Type1_doy2020001_aid0001
name : MCD12Q1.061_LC_Type1_doy2020001_aid0001
min value : 1
max value : 17
raster): anthropogenicAgain using the package geodata.
class : SpatRaster
dimensions : 2400, 1680, 1 (nrow, ncol, nlyr)
resolution : 0.008333333, 0.008333333 (x, y)
extent : -82, -68, -19, 1 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source(s) : memory
varname : wildareas-v3-2009-human-footprint_geo
name : wildareas-v3-2009-human-footprint
min value : 0
max value : 46
Rasters should have the same spatial extent and resolution… If they don’t match… what should we do?
xmin ymin xmax ymax
-81.5 -18.5 -68.5 0.0
xmin ymin xmax ymax
-81.5 -18.5 -68.5 0.0
xmin ymin xmax ymax
-81.5 -18.5 -68.5 0.1
xmin ymin xmax ymax
-82 -19 -68 1
xmin ymin xmax ymax
-82 -19 -68 1
xmin ymin xmax ymax
-81.4109426 -18.3479754 -68.6650797 -0.0572055
We will use the package terra to process the rasters.
First we need to create a new raster to use as a template.
Bug in terra
We could do:
r <- rast(res = 50000,
ext = extent, # or even `ext(vect(peru))`
crs = "epsg:4326")
But there’s a bug :O
I reported it https://github.com/rspatial/terra/issues/1618
We will use the package terra to process the rasters.
First we need to create a new raster to use as a template.
r <- rast(res = 50000, # Target cell size
ext(vect(st_transform(peru, proj_peru))), # Extent of the raster
crs = proj_peru) %>% # Useful planar projection
project(., "epsg:4326") # Projection of our rasters
rclass : SpatRaster
dimensions : 41, 29, 1 (nrow, ncol, nlyr)
resolution : 0.4514065, 0.4514065 (x, y)
extent : -81.74007, -68.64928, -18.49128, 0.01638293 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
We will resample() all of the rasters to our target resolution.
And finally we will use mask() to intersect the layer with Peru.
We will use the package tmap.
Annual mean temperature (WorldClim BIO1)
Precipitation seasonality (WorldClim BIO15)
Elevation (WorldClim)
Land-cover (MODIS Terra MCD12Q1). We need to make sure values are factors.
Then we can add values to the levels.
levels(land[[1]]) <- data.frame(
ID = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17),
label = c(
"Evergreen Needleleaf Forest",
"Evergreen Broadleaf Forests",
"Deciduous Needleleaf Forests",
"Deciduous Broadleaf Forests",
"Mixed Forests",
"Closed Shrublands",
"Open Shrublands",
"Woody Savannas",
"Savannas",
"Grassslands",
"Permanent Wetlands",
"Croplands",
"Urban and Built-up Lands",
"Cropland/Natural Vegetation Mosaics",
"Permanent Snow and Ice",
"Barren",
"Water Bodies"
)
)You might not have the same classes/levels at different areas. Check here the classification values: MCD12_User_Guide_V6 (see page 7 Classification Legends).
Land-cover (MODIS Terra MCD12Q1)
Human footprint index (SEDAC)
Let’s get values per grid-cell.
[1] "PER_wc2.1_30s_tavg_1" "PER_wc2.1_30s_tavg_2" "PER_wc2.1_30s_tavg_3"
[4] "PER_wc2.1_30s_tavg_4" "PER_wc2.1_30s_tavg_5" "PER_wc2.1_30s_tavg_6"
[7] "PER_wc2.1_30s_tavg_7" "PER_wc2.1_30s_tavg_8" "PER_wc2.1_30s_tavg_9"
[10] "PER_wc2.1_30s_tavg_10" "PER_wc2.1_30s_tavg_11" "PER_wc2.1_30s_tavg_12"
Let’s see all the values together.
Let’s see all the values together.
tavg prec elev land
Min. : 2.596 Min. : 1.62 Min. : 15.58 Min. : 1.000
1st Qu.:14.045 1st Qu.: 807.77 1st Qu.: 196.75 1st Qu.: 2.000
Median :23.439 Median :1710.85 Median : 660.72 Median : 2.000
Mean :19.666 Mean :1625.41 Mean :1483.44 Mean : 6.414
3rd Qu.:25.921 3rd Qu.:2418.22 3rd Qu.:2779.24 3rd Qu.:10.000
Max. :26.983 Max. :4495.59 Max. :4686.61 Max. :17.000
NA's :612 NA's :612 NA's :661 NA's :573
footp
Min. : 0.000
1st Qu.: 0.000
Median : 2.000
Mean : 3.002
3rd Qu.: 4.000
Max. :34.000
NA's :573
Let’s see how these values correlate.
What can this plot tell you?
Let’s read the data, plot it and see the patterns
Let’s read the data, plot it and see the patterns
Species richness vs mean annual temperature
Species richness vs precipitation seasonality
Species richness vs elevation
Species richness vs land-use
Species richness vs human footprint
raster)
We are done! No it’s your turn :)
raster)This time, you will use the following drivers:
See the solutions at practical_week3_STUDENTS_FILE_solved.qmd.